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Summary. A new code for evaluation of light absorption and scattering by in- 
terstellar dust grains is presented. The radiative transfer problem is solved using 
ray-tracing algorithm in a self-consistent and highly efficient way. The code demon- 
strates performance and accuracy similar or better than that of previously published 
results, achieved using Monte-Carlo methods, with accuracy better than ~ 3% for 
most cases. The intended application of the code is spectrophotometric modelling 
of disk galaxies, however, it can be easily adapted to other cases that require a 
detailed spatial evaluation of scattering, such as circumstellar disks and shells con- 
taining both point and distributed light sources. 



1 Problem statement 

The purpose for the developing radiative transfer problem solving code, de- 
scribed in this article, was to model spatial and spectral energy distribution 
(SED) observed in external galaxies. The nature of this problem requires 'self- 
consistency' of a solution - the resulting SED of a model must depend only 
on the SED of stellar sources and assumed properties of dust without any 
preconditioning on light and attenuation distribution within galaxy [TVA03] . 

While galaxies in general are complex objects with three-dimensional (3D) 
distribution of radiation and mater, in most cases they are dominated by ax- 
ial symmetry (2D), allowing significant simplification of the model geometry. 
However, the model should account for presence of macroscopic structure 
within galaxies, possibly including elements having other symmetry, such as 
bars and spiral arms (2D+). 

Most present day astrophysical radiative transfer codes employ either a 
Monte-Carlo (MC, eg. [CF01]) or a ray-tracing (RTR, eg. [RS99], [RM02]) 
methods. Some of implementations of these methods were compared by [BD01] 
for ID and by [DT02] for 2D cases. The RTR approach allows the optimiza- 
tion of solution for a given system geometry, which was the main reason to 
use it as a basis for the Galactic Fog Engine (hereafter 'GFE'), a program for 
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self-consistent solution of radiative transfer problem in dusty media with pri- 
marily axisymmetric geometry. This paper concentrates on radiative transfer 
in ultraviolet-to-optical wavelength range assuming exclusively coherent scat- 
tering, therefore in most equations the wavelength dependence will be omitted. 

2 Algorithm 

2.1 Model description 

The GFE iteratively solves a discrete bidirectional radiative transfer problem, 
producing intensity maps of the model under arbitrary inclination at a given 
wavelength set. The foundation of the iterative evaluation of radiation transfer 
equation was laid out by Henyey [Hen37]. By solving a set of one-dimensional 
radiative transfer equations 



where <P and u denote the scattering phase function and albedo and k and 
j are absorption and emissivity coefficients of the medium, the initial system 
SED is separated into: escaped energy, that reaches an external observer; 
energy, absorbed by grains, that is eventually emitted as thermal radiation; 
and scattered energy distribution. The solution is then repeated, substituting 
scattered energy as initial distribution for the next iteration, accumulating 
resulting escaped and absorbed energy, until certain convergence criteria are 
met, either a fixed number of iterations, or remaining scattered energy being 
below specified threshold. After convergence is reached the dust temperature 
is calculated from absorbed energy distribution. If it is necessary to account 
for self-scattering of thermal radiation by dust grains, the resulting emission 
SED can be input into scattering evaluation loop and the process repeated 
until the final radiative energy distribution is obtained, and then used to 
produce SED as seen by an external observer. 

Calculations are performed within a cylinder with a radius r m and height 
above midplane z m , which is subdivided into a set of layers of concentric, 
internally homogeneous rings ('bins') of arbitrary radial and vertical thick- 
ness. Since the linear extent of each individual light source (star) is negligible 
compared to the size of system, it is possible to solve the radiative transfer 
problem considering every volume element of the model having both attenu- 
ating (light scattering and absorption by interstellar dust) and emitting (light 
emission by the stars and thermal radiation of the dust particles) properties 
per unit volume, defining for each ring denoted by indexes r and z its total 
absorption k' — n(r, z), and emissivity j' — j(r, z, a, S), combined from inter- 
nal light sources and energy scattered within its volume, with angles a and 5 
defining the direction of radiation propagation. 
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Fig. 1. The model geometry. Panel a) shows the diametral, panel b) - central 
plane cross-section of the model. The distribution and density of ray-tracing paths 
(shown as arrows) are computed to produce even and sufficient sampling of the model 
volume. Panel c) illustrates the discrete radiative transfer in cross-section parallel 
to the model Z axis. Limits of plane-parallel layers for one-directional treatment 
are shown as dotted lines while boundaries of the individual rings are outlined with 
solid lines. 



2.2 Radiative transfer equation 

GFE uses static ray-casting geometry, determining the set of rays that ensures 
a required degree of sampling of the model volume (fig. la and lb). If the 
viewing solid angle, containing each direction, can be held small, the radiative 
transfer along these rays can be solved as in plane-parallel homogeneous layer 
case for a series of intervals traversing rings until crossing the outer boundary 
of the model (fig. lc). In a form suitable for computer implementation an 
incident intensity on a given point for a light path separated into n intervals 
of length k, numbered outwards from that point, is written as 

^-±(ne- k 'A§(l^e-^). (2) 
i=i \j=i J 2 

Similarly, the intensity of radiation scattered with albedo ui from a given 
direction (a : 5) into all other directions within a certain interval denoted by 
index "1" is 
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(3) 



When considering azimuthally inhomogeneous model configuration (3D 
case), each ring is subdivided into required number of azimuthal segments. 
The number and directions of rays cast through the system have to be modified 
accordingly to include new sets of rays in azimuthal direction, however, the 
ray-tracing part of the algorithm is unchanged. The computational time scales 
as x logiVbin for 2D and A^/ 3 x log7V bi n for 3D cases. 

2.3 Scattering phase function 

Since angular distribution of radiation at each point in the model is non- 
isotropic, it must be described using a numerical phase function (matrix), 
providing the radiation intensity towards a set of predefined reference di- 
rections ('RDs') described by angular coordinates (ao : So)- There exists a 
number of ways to distribute RDs on a sphere, however, those methods that 
produce a set of RDs arranged in iso-latitude rows are the most efficient in 
this particular model geometry, allowing both efficient storage and retrieval of 
scattered intensity and fast rotation of the phase matrix around model Z-axis. 
The memory requirements and the overall algorithm's performance have also 
to be taken into consideration. 

In this work the following methods of RD distribution were compared: 
HEALPix 1 [GHW99], HTM 2 [KST01], a trivial iso-latitude triangulation 
(hereafter 'TT'), fig. 2a) and a square matrix with elements ('texels') cor- 
responding to evenly spaced (ao : So) coordinates (hereafter 'Texel'). For 
triangulation schemes and HEALPix the radiation intensity towards a given 
point was interpolated between 3 nearest RDs using either 'flat' (fig. 2b) or 
'spherical' (fig. 2c) weights. 

3 Computational precision 

3.1 Scattering phase function interpolation 

To compare used sphere subdivision and interpolation algorithms a following 
test model (hereafter a 'standard model') was employed: a cylinder with height 
to radius ratio z m /r m = 0.2, divided into A^ bin = 441 (21 x 21) equally spaced 
rings, filled with radiating and absorbing particles whose density follows a 
double exponential law 



x http : / /www . eso . org/science/healpix/ 
2 http: //www. sdss . jhu.edu/htm/ 



Radiative 



transfer problem in dusty galaxies 5 



a) 




Fig. 2. The reference point structure used for interpolation of the scattering phase 
function. Panel a) shows the trivial iso-latitude triangulation for one hemisphere 
with reference directions arranged symmetrically against diametral and horizontal 
planes. Panels b) and c) represent two implemented interpolation schemes, 'flat' 
and 'spherical'. In case of spherical interpolation, input from each triangle vertex is 
weighted by the area defined by shortest distances from the given direction to the 
vertices (Si for 1-st vertex and so on). 



p(r,z) = poe- r/r "e- z ^ zo (4) 

with r = 0.2r m and z n — 0.2z m and the model central optical depth per- 
pendicularly to the central plane T ct = 25. The Henyey & Greenstein [HG41] 
scattering phase function 

1 - Q 2 
(1 + g 2 -2gcos9) 3/2 

was used with asymmetry parameter g = 0.75. 

The primary quality criteria of a given algorithm is its ability to represent 
the angular intensity distribution of anisotropic scattering. If the phase func- 
tion representation is exact, the distribution of values (<P'(6) — <P(9))/<P(9) 
(where <P'(8) is a resulting numerical phase matrix) would be a 8- function. 
However, since employed methods introduce different types of numerical er- 
rors, the actual distribution form depends strongly on the phase function 
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sampling and the interpolation algorithm. An examples of resulting error dis- 
tributions (as relative numerical phase matrix deviation from its analytical 
form) for the algorithms tested are shown on fig. 3. 

As can be seen, methods providing uniform sphere coverage produce more 
preferable error distributions. With the increasing number of RDs the repre- 
sentation of the scattering phase function improves, reducing maximal possi- 
ble deviation from the true value, particularly for TT (fig. 3a) and HEALPix 
(fig. 3c) methods, with TT algorithm showing slightly better error distri- 
bution form. 'Spherical' interpolation scheme (fig. 3, right column in panels 
a - c) produces more symmetrical error distributions, while 'flat' interpola- 
tion (left column) in some cases introduces additional numerical errors. When 
compared with other methods, attempt to reproduce scattering phase function 
using simple matrix with no interpolation between its elements (Texel scheme, 
fig. 3d) produces results of average quality, its error distribution quickly reach- 
ing a 'saturated' form with increasing number of RDs. The described error 
distribution is somewhat dependent on the orientation of the scattering phase 
function relative to the set of RDs, this dependence being minimal for the 
methods with identical size of interpolation elements (HEALPix). All meth- 
ods that use interpolation display similar performance for a given number of 
RDs (table 1). 

Table 1. The normalized computing time for models using different sphere subdi- 
vision and scattering phase function interpolation algorithms. 



Interpolation 


HEALPix 


HTM 


TT 


Texel 


method 










'flat' 


2.5 


2.8 


2.5 


1.0 


'spherical' 


8 


10 


8 





3.2 Volume sampling and subdivision 

The problem encountered applying numerical methods is error accumulation. 
In case of iterative ray-tracing it arises from sampling and interpolation er- 
rors. The source of sampling errors is incomplete/inadequate spatial sampling 
of the system while interpolation errors are related to the scattering phase 
function approximation, described in the previous section. Both error types 
independently affect every ray traced through the system, thus the accumu- 
lated error increases with the increasing number of bins and rays. This makes 
oversampling undesirable not only due to increasing computational time, but 
also for a reason of minimizing numerical errors. 

As a measure of method's quality a defect in energy balance E CII as a 
percentage of total energy radiated within system E tot 
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Fig. 3. The distribution of relative numerical phase matrix deviation from its an- 
alytical form for different sphere subdivision algorithms. Panels a) - d) correspond 
to TT, HTM, HEALPix and Texel methods. For the first three methods the results 
obtained using both 'flat' (left panels) and 'spherical' (right panels) weighted inter- 
polation are presented. Thin line shows the results obtained for approximately 3100, 
thick line - for approximately 12000 reference directions. 
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t~, -E'tot -E a bs -E S ca ^csc fn \ 
#crr = = (6) 

^tot 

is used. Here E esc , E a b s and E sca are the parts of a total radiated energy that 
escaped the system, was absorbed and remained to be scattered within the 
system, respectively. 

To determine the sampling and gridding influence on the model precision 
the following two tests were performed. Firstly, the radiation field in the stan- 
dard model using TT algorithm with scattering phase function represented 
by a set of 182 RDs was computed with different number of rays N ra y , cast 
through each ring. The results, presented on fig. 4a, show a significant error 
accumulation effect. Then, keeping a number of rays per ring constant the 
number of rings iVbi n in model was changed (fig. 4b). As can be seen, improv- 
ing the sampling of a model decreases the approximation errors to a certain 
minimum, limited by internal errors of a chosen scattering phase function in- 
terpolation method. This geometric configuration can be considered optimal, 
since with further increase in a number of bins and rays the quality of the 
solution begins to deteriorate due to error accumulation. 
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Fig. 4. The influence of subdivision and sampling of the model on the energy 
balance. Panel a) displays the energy defect for a given number of rays cast per 
model bin; panel b) shows the same quantity for a models consisting of different 
number of bins. 



3.3 Dust optical properties 

Other important aspect of a numerical radiative transfer solution is its sen- 
sitivity to variations in scattering parameters: albedo u and asymmetry pa- 
rameter g. Model precision and stability for different w and g values place 
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a constraint on the wavelength range where a given method can be applied. 
The influence of scattering parameters on energy defect was analyzed using 
the standard model with 7V bin = 441 and r ct = 10. 




Fig. 5. The dependence of the energy losses within model on the grain scattering 
parameters: scattering phase function asymmetry g (panel a), with albedo assumed 
to be lu — 0.5 for all g values), and albedo lj for g — 0.6 (panel b). 



The dependence of overall model precision on scattering asymmetry pa- 
rameter g is shown on fig. 5a. The total energy defect after 9 iterations show 
some variation with the increasing g up to the limit imposed by the angular 
scattering phase function gridding (182 RDs) used in the calculations, after 
which the energy losses in the scattering phase matrix render results invalid. 

The effects of the grain albedo on the model accuracy and stability were 
investigated using similar method. All computations were performed for 7 
iterations assuming g = 0.6. The results are shown in fig. 5b. Within the 
range of lj values, applicable to astrophysical dust grains, those errors stay in 
acceptable limits, and do not influence the stability of the solution. 



4 Summary 

The code described in this paper has undergone an extensive testing and shows 
the flexibility and performance satisfying the requirements for the models of 
the global radiation transfer in dusty galaxies [SV02]. It has been success- 
fully applied to model both integral and position dependent SEDs of several 
galaxies, some of the first results presented in [SSV03]. 

The main limiting factor affecting the applicability of the described code 
is the scattering asymmetry parameter g\. In order to correctly treat the 
scattering with g\ approaching 1, the number of required reference directions 
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rises sharply, affecting both the performance and the precision of the method. 
For a disk galaxy model such as one decribed in this paper a satisfactory 
convergence is obtained for g\ in [ -0.8 ; 0.8 ] range which includes the optical 
properties of typical astrophysical grains scattering photons from microwave 
up to extreme UV wavelengths. 

Other model properties, such as optical depth t\ and the relative amount 
of scattered radiation (dependent on albedo u\) seem to have relatively lit- 
tle effect on the quality of the solution. However, models with large optical 
depth (of order of a few 100's), particularly having a steep matter distribution 
gradients, require a significant amount of computing time and storage. 

The application of this code is not restricted to the systems with dis- 
persed sources and absorbers, the algorithm being easily extended to include 
treatment of interaction between radiation field and surfaces of macroscopic 
objects. 
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